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H ' Abstract 

lyj ■ In intensity-modulated radiation therapy, optimal intensity distribu- 

O ' tions of incoming beams are decomposed into linear combinations of leaf 

' jy^ , openings of a multileaf collimator (segments). In order to avoid inefficient 

^~j' dose delivery, the decomposition should satisfy a number of dosimetric 

,S^ [ constraints due to suboptimal dose characteristics of small segments. How- 

C^< ever, exact decomposition with dosimetric constraints is only in limited 

cases possible. The present work introduces new heuristic segmentation 

fvj , algorithms for the following optimization problem: Find a segmentation 

^1 of an approximated matrix using only allowed fields and minimize the 

QQ ' approximation error. Finally, the decomposition algorithms were imple- 

^\ , mented into an optimization programme in order to examine the assump- 

OO ' tions of the algorithms for a clinical example. As a result, identical dose 

^^ ' distributions with much fewer segments and a significantly smaller num- 

|/-\ . ber of monitor units could be achieved using dosimetric constraints. Con- 

f*"^ ' sequently, the dose delivery is more efficient and less time consuming. 
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1 Introduction 

In intensity-modulated radiation therapy (IMRT), intensity matrices with non- 
negative integer entries are computed for each irradiation field. After dis- 
cretization of the field into bixels, each entry of the matrix corresponds to the 
required intensity within this bixel. The segmentation step consists in decom- 
posing the matrix into a linear combination of subfields (segments) shaped by 
a multileaf collimator (MLC). The first intuition is that a treatment plan is opti- 
mal, if the linear combination of the chosen segments equals the matrix. Such 
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Figure 1: Electron dose output at the dose maximum normalized to the dose 
output of the 10 cm x 10 cm field and electron penetration depth of the 90 % 
depth-dose as a function of square field size and electron energy (from fTS]). 
The fields were shaped by an add-on MLC for electrons presented in Figure 
|6l A minimum MLC field size of approximately 3 cm x 3 cm is necessary 
for decomposing intensity distributions into leaf openings to ensure an output 
factor of nearly 1 and an energy-dependent penetration depth. 



a plan consists of various segments possibly including those segments where 
most of the irradiation field is covered and only few bixels receive radiation. 

For dosimetric reasons, however, the model assumption is not given in 
practice. Irradiation of small photon or electron segments result in a much 
lower dose output compared to conventional conformal fields. Therefore, the 
linearity assumption, that irradiating one segment is equivalent to dividing it 
into two parts and irradiating them separately, only holds, if the two parts are 
still sufficiently large. In addition, the penetration depth of electrons decreases 
with decreasing field size and is almost independent of the beam energy for ap- 
proximately 1 cm X 1 cm fields. However, the energy dependence of the pene- 
tration depth is necessary for our new IMRT technique with electron beams to 
adjust the dose to the target volume by use of various beam energies. Figure [1] 
shows that electron fields of approximately 3 cm x 3 cm are necessary to keep 
an output factor of nearly 1 and an energy-dependent penetration depth. 

As a consequence, a treatment plan should consist of segment shapes satis- 
fying certain constraints that ensure a minimum field size. For practical pur- 
poses it is also necessary that the field openings are connected and do not de- 
generate into two or more parts. Besides those dosimetric constraints, there 
are also technical constraints reducing the number of allowed shapes. One is 
the leaf overtravel constraint that accommodates the fact that the left (respec- 
tively right) leaf of the MLC cannot be shifted further than a threshold to the 
right (respectively left). These constraints have the consequence, that not every 
intensity matrix is decomposable in segments satisfying the constraints. This 
leads us to the task to find an approximation matrix and its decomposition into 
"good" segments, that differs from the given intensity matrix as little as pos- 
sible. The aim is to generate equivalent treatment plans with good segments 
leading to a reduction in the segment number and monitor units, respectively. 

The decomposition problem for the exact case without concerning any ad- 
ditional constraints is well studied. Algorithms for the minimization of the 
beam-on time can be found in l2ll3ll6l [T4l[T5ll2TI . Approaches for minimiz- 
ing the number of used segments are given in ll7l [T8ll26l . A variety of techni- 
cal constraints are considered, see fS" 161 for the interleaf collision constraint, 
that prohibits an overlap of adjacent leaf pairs, and L16l [TTl l22l l23l [27.1 for the 



tongue-and-groove constraint. Kamath et. al. Il2l1 also investigate the mini- 
mum separation constraint that requires a minimum leaf opening in each row 
and develop a criterion for a matrix being decomposable under this constraint. 
Engelbeen and Fiorini flOl deal with the interleaf distance constraint where 
the allowed difference between two left (respectively right) leaf positions is 
bounded by some given threshold. 

An approximation problem with the aim of reducing the total beam-on 
time was first formulated in |8| and generalized to approximated decompo- 
sition with interleaf-collision constraint in [19J and |20|. The dependence be- 
tween field size and output factors, penetration depth and depth-dose fall-off 
is outlined in fT3). These considerations lead to the decomposition problem 
using segments that satisfy some minimum field size constraints. Under these 
constraints, an exact decomposition of the intensity matrices is, in general, no 
longer possible (cf. [21]) and an approximation problem has to be formulated. 

Another algorithmic approach that aims at minimizing the number of seg- 
ments while keeping the quality of the treatment plan is the direct aperture 
optimization that combines the choice of beams, apertures and weights with- 
out computing a leaf sequencing step. Shepard et al. |28| allow only a limited 
number of apertures for each beam, Bedford and Webb [4J also integrate con- 
straints on the segment shape and size in the direct aperture optimization ap- 
proach. Our algorithm is applicable if one uses intensity profile segmentation 
and wants to compute segmentations satisfying certain field size constraints 
and reducing the complexity of the plan. Matuszak et al. [25 1 deal with the 
minimization of the monitor units by smoothing the intensity profiles. 

The paper is organized as follows. Section |2] gives two definitions of what 
we call a segment with good dosimetric properties, one basic definition and 
an extended one including one further constraint. We concretely define the 
approximation problem and in Section|3]propose heuristic algorithms for both 
definitions, each of them consisting of seven different steps. The different parts 
of the algorithm and their properties are analyzed in Section |4l We especially 
outline, that the solutions of the subproblems in step 1 and 2 are indeed opti- 
mal. Section|5]introduces the clinical case we used for testing the quality of our 
segmentations. Section |6] gives computational results for the test case and de- 
tailed numerical results for the segmentation of clinical matrices from different 
IMRT treatment plans. 

2 Problem formulation and definitions 

Throughout the paper we use the notation 

[k] = {l,2,...,k} and [k,l] = {k,k + 1,. . . ,1} 

for integers k and I, k < I. Let A = (a^) denote the given fluence matrix of 
size m X n. Feasible leaf positions of the MLC are modeled as binary matrices 
S = (sij), called segments, that satisfy the consecutive-ones-property in each 
row. In other words, S" is a segment, if there are integral intervals [ii, r^] for all 
i e [m], representing the positions of the left and the right leaf, such that 
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1 if k <j< n 



((*,j) e [m] X [n]). (1) 



I otherwise 
Furthermore, for each segment S, we define sqj — s.,„+i ^ = for all j G [n\. 



For the described reasons, we introduce five parameters bi, br, gi, §2 and / 
representing the following constraints: 

(i) Left Leaf Overtravel Constraint: For all i G [m], we require k < bi. In 
each row the left leaf cannot be shifted more to the right than to the bixel 
with index bi . 

(ii) Right Leaf Overtravel Constraint: For all i e [m], we require r^ > br- In 
each row the right leaf cannot be shifted more to the left than to the bixel 
with index br- 
im) Minimum Separation Constraint and Row Overlap: If a row i G [m] is 
not totally covered, we require Vi — U > 51 — 1. Similarly, if rows i and 
i + 1 are not completely covered, we claim min(ri, r^+i) — niax(Zi, U+i) > 
gi — 1. At least gi consecutive bixels in each row receive radiation and 
the irradiated area of two consecutive rows overlaps in at least gi bixels. 

(iv) Minimum Vertical Gap: We require a minimum vertical field size and a 
minimum vertical size of the covered regions, i.e. in each column con- 
secutive ones or zeros should have a minimum number 92- In detail, if 
Si^ij — 0, Sij — Si+i.j = • ■ • = Sk-ij = 1 and Skj = for some column 
j, we have k — i > g2. Analogously, we require the same for consecutive 
zeros framed by ones. 

(v) Minimum Total Field Height: At least / consecutive rows of the field are 
not totally covered, i.e. there are at least / consecutive rows with / < r. 
This ensures, that the total size of the field is reasonably large. 

Of course, these parameters only make sense iil < bi,br,gi < n,l < g2, f < m 
as well as 1 < br < bi < n and 52 < ./• 

The case that one row of the field is totally covered and receives no radiation 
at all, is throughout this paper represented by the leaf positions I = n + 1 and 
r = 0. In practice, one will of course choose leaf positions of the form I = r + 1 
with I < bi and r > br that respect the leaf overtravel constraints (i) and (ii). 

Remark 1. The Minimum Vertical Gap can be formulated in terms of the 
leaf position as follows: If k < k-i for some i > 2, then we also require k+i < 
h,li+2 < h, ■ ■ ■ ,li+g2-i < h- Analogously, if Ti > r^-i for some i > 2, we 
also have r^+i > r^, ri_|_2 > ^i, • • • , '^1+92-1 ^ ^i- This ensures, that in vertical 
direction, we always have at least 172 bixels open. Therefore, ^ve additionally 
not allow that Ii > k^i or r^ < r^-i for 2 < i < g2 and forbid also k < k^i 
as well as r^ > r^-i for m — g2 + l<i<Tn. Similarly, we require at least 
52 bixels closed in vertical direction, if there are open bixels above and below 
in this column of the matrix. Thus, we make sure that also thin shapes in 
vertical direction, having negative dosimetric properties as discussed in the 
introduction, are forbidden. 

A segment is called connected if the irradiated area that corresponds to its 
leaf positions does not resolve into two or more parts, i.e. if the corresponding 
rectilinear polygon (considered as an open set) is connected. 

As the realization of the minimum vertical gap turns out to be the most 
difficult task, we introduce two different definitions of "good" segments. 

Definition 1. Given the parameters fe(,&r,3i,/ G N with 1 < bi,br,gi < n, 
1 < / < m as well as I < br < bi < n, a segment with good dosimetric properties S 
is a connected segment satisfying the constraints (i), (ii),(iii) and (v). 



Figure 2: For gi = 52 = 2 the two left segments violate the minimum vertical 
gap constraint and only satisfy (iii), whereas the two right segments meet both 
conditions (iii) and (iv). 



Definition 2. Given the parameters bi,br,gi,g2, f G N with 1 < bi,br,gi < n, 
1 < 52 , / < m as well as 1 < br < bi < n and g2 < /, a segment with very good 
dosimetric properties is a connected segment satisfying the constraints (i)-(v). 

For brevity of notation we will call the segments with good dosimetric 
properties from now on simply segments and the segments with very good 
dosimetric properties advanced segments. 

All in all, we have two optimization problems: Given a matrix A with posi- 
tive integer entries and suitable parameters bi,br,gi,g2, f, find a segmentation 



B = y^^UkSk 



fc=i 
where the Uk are positive integers and 

• Approximated Segmentation into segments (ASS): the Sk are segments 

• Approximated Segmentation into advanced segments (ASAS): the Sk 

are advanced segments 



such that 






The value of the objective function of the optimization problem is called 
total change. The delivery time of the segmentation is X)fc=i ^fe ^^^^ the number 
of segments is t. 

Obviously, for the parameter set bi — n, br — 1, 51 = .92 = / = 1/ segments 
are simply connected segments in the sense of Equation (|1), the approximation 
problem has as value of the objective function and degenerates to a segmen- 
tation problem into connected segments defined by Equation |(T]|. 



3 Approximated segmentation 

Now we introduce two basic algorithms for ASS and ASAS, each consisting 
of seven different steps. Two of the steps are identical in both algorithms as 
they are computed before the segmentation step and therefore do not affect 
the parameter 92- Five of the steps differ subject to whether the constraint 92 



is regarded or not (namely steps 3-7). Steps 1 and 2 will be solved exactly, 
whereas steps 3-7 are heuristic. 

The basic structure of our algorithms is as follows: 

1 . Solve the Leaf Overtravel Constraint Problem (LOC) on A: Given bi,br G 
N with I < bi < br < n, find an approximation B with nonnegative in- 
teger entries that can be decomposed with respect to the Leaf Overtravel 
Constraint such that the total change is minimal. 

A =8. 

2. Solve the Minimum Separation Constraint Problem (MSC) on A: Given 
5i G [n], find an approximation B that can be decomposed into segments 
with Ti — li> gi — 1 for all i G [m] such that the total change is minimal. 

A:=B. 

3. ASS: Compute an approximated segmentation S into connected segments 
satisfying (i)-(iii), but not necessarily (v). 

ASAS: Compute an approximated segmentation S into connected seg- 
ments satisfying (i)-(iv), but not necessarily (v). 

The approximated segmentation may violate (v) and belongs to an ap- 
proximation matrix B that might have a large total change. 

4. Combination of Fields: Combine stepwise two disjoint fields, if this is 
possible with small total change. For ASAS, make sure that the new field 
satisfies (iv). 

5. Make-two-of-two: For S, S' G S, where S violates (v) and S" satisfies (v), 
compute a substitution 5 + S" = S + S' , such that S and S' satisfy (v). For 
ASS, make sure that S and 5" satisfy (i)-(iii). For ASAS, S and S' must 
also satisfy (iv). 

6. Handle Critical Segments: If there are still segments violating (v), try to 
combine them with a feasible segment such that the total change of the 
combination is smaller than omitting the critical segment. Take care, that 
(i)-(iii) (resp. (i)-(iv)) hold. 

7. Total Change Improvement: For all segments and all rows i, check whether 
an increase or decrease of U (respectively r^) reduces the total change. For 
ASS, change U or r^ only if (i)-(iii) still hold. For ASAS, change U or r^ 
only if (i)-(iv) still hold. Look at all the segments cyclically until no more 
changes are possible. 

The output of the algorithm is a segmentation of an approximation ma- 
trix that consists only of segments for ASS and only of advanced segments 
for ASAS. The LOC- and the MSC-approximation aim at producing a first ap- 
proximation that can be better decomposed with the given constraints (i)-(iii) 
in the segmentation step than the initial matrix. The combination of fields, the 
make-two-of-two-step and the handling of critical segments try to provide seg- 
ments satisfying (v) without producing too much total change. Finally, the total 
change improvement is computed in order to improve the approximation. 

A very important feature of our heuristic is that further constraints can eas- 
ily be taken into consideration. For example, if adjacent left and right leafs may 
not overlap (interleaf collision constraint) one can allow only leaf positions that 
respect this constraint in the optimization steps 3-7. 



4 Subproblems 

Now we describe in detail the subalgorithms and outline some of their basic 
properties. Throughout the steps that follow the leaf overtravel approximation, 
we allow only leaf positions that respect the leaf overtravel constraint. For 
simplicity, we will not mention this basic fact in every step. 

4.1 Leaf Overtravel Constraint 

As the leaf overtravel constraint only affects a single row of the matrix, the 
problem LOC can be solved for each row independently Thus, we compute an 
optimal approximation of a vector a. Segmentations reduce to sums of inter- 
vals [l,r]. Segments are simply 0-1 -vectors s with consecutive ones. 

k 

Lemma 1. A vector a has a segmentation a = J2 ^i ^'^^ corresponding leaf positions 

1=1 
li < biandri > briffaj > aj+iforall j G [bi,n — l]andaj > aj-iforallj G [2,6^]- 

Proof. Let gq = a„+i := and a+ := max{0, a}. On the one hand, the algorithm 
of Bortfeld (see 161) provides a segmentation where the left leaf position is j 
for exactly {gj — aj_i)+ segments. Analogously, the right leaf position is j for 
(flj — flj+i )+ segments and no other leaf positions occur. On the other hand, it is 
obvious that if gj > gj^i (respectively gj > Gj+i) there will be a segment with 
left (respectively right) leaf position j in every segmentation. This concludes 
the proof. D 

Therefore, we have to find an approximation vector, that has no up-steps 
after index bi and no down-steps before index br- As we assume b,. < bi, we can 
use symmetry to solve the approximation problem for the right leaf positions. 
Besides, the criterion from Lemma [l] shows, that bj = gj for j G [br + l,bi — l] for 
each optimal solution of LOC. We simply need to solve the following problem 
for the subvector (afc, , . . . , a„) and the left overtravel constraint: 

LOC-left: Given a vector v = {vi, . . . ,Vk), find an approximation vector w 

k 

with Wj > u'j+i for j G [A: — 1] such that || ?; — to || i= ^ |wj — Wj | — > min. 

The algorithm for solving the problem LOC-left is described in Algorithm[T] 
in the appendix. It uses a graph theoretical approach and computes a shortest 
path in a layered digraph, where the j-th layer consists of nodes representing 
the possible entries of the j-th component of the approximation vector. The 
problem LOC-left is similar to the Monotone Discrete Approximation Problem 
(MDAP) formulated in |8| and the algorithm follows the same idea. 

Let min := min^grj.] Vj and mGX := maxjgr^ji Vj and let tcij be the objective 
value of an optimal solution of LOC-left with Wj = i. Let prejj be the corre- 
sponding predecessor Wj - 1 . With respect to Algorithm[T](that uses the notation 
above) we yield the following 

Theorem 1. Algorithm\l\computes an optimal solution of LOC-left. 

Proof. The initial values tcii are trivially correct. Let now j > 1 andlet {wi, ... ,Wj) 
be an optimal approximation of (ui, ... ,Wj) with Wj = i. By induction, icm„jj_i 
is computed correctly and thus 

i 

^\VI -Wi\= tCnij_i,j-l + \Vj -i\> tCij. 

1=1 



Therefore tcij is a lower bound for the total change. The choice of iopt makes 
sure, that the optimal value of Wk is chosen and obviously, the approximation 
vector from Algorithm[T]realizes the lower bound for the total change of tci^^^ k ■ 

op .^ 

4.2 Minimum separation constraint 

Like the Overtravel Constraint, the Minimum Separation Constraint can be 
handled independently for each row of the matrix. Thus the task is the foUow- 



MSC-Row Given a vector v ~ (wi , . . . , u„) with nonnegative integer entries, 
find an approximation vector w with nonnegative integer entries, such that w 

n 

has a decomposition into intervals of length > gi and || ?; — u? ||i= X) l^j ~ 

Wj I -^ min. 

We know from Kamath et al. (see [21]) that a vector a can be decomposed 
without violating the minimum separation constraint, if the optimal decompo- 
sition of their algorithm SINGLEPAIR does not violate the minimum separa- 
tion constraint. For example, the vector a — (1, 2, 1) cannot be decomposed 
with gi = 3, as the optimal decomposition is (1, 2, 1) = (1, 1, 0) + (0, 1, 1) and 
the used intervals do not have a minimum length of 3. This motivates the ap- 
proximation problem defined above. 

Obviously, the problem MSC-Row can be formulated as an integer linear 
programming problem as follows: 
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We solve this integer program for each row of the matrix using SCIP |T| with 
SoPlex 1 29 1 as LP solver. The problem can also be solved by a combinatorial 
algorithm using a minimum cost flow formulation which is shown in LllJ . 

4.3 Segmentation 

Let B be the approximation matrix resulting from the MSC-step. The basic 
idea of the segmentation is to consider the current total change to be the sum 
of the absolute values of the entries of B and to iteratively compute a segment 
S whose subtraction reduces the current total change. In each step, the matrix 
B is updated by setting B := B — S . At the end of the segmentation, a positive 
entry in B represents a bixel with underdose and a negative entry a bixel with 
overdose. 



4.3.1 Segmentation for ASS 

Let bi denote the i-th row of B for all i G [m]. A segmentation consists of seg- 
ments S each represented by its leaf positions li, I'i for i G [m]. The main body 
for the segmentation step is described in Algorithm |2] in the appendix. This 
algorithm uses the subroutine Find interval ASS that is precisely described in 
Algorithm |3l 

The idea behind this heuristic choice of the segment S being subtracted 
from B in each step is, that we compute the first interval from a sliding window 
segmentation (see again [6]), with / as index of the first up-step and r as index 
of the first down-step in the corresponding row. If these values already satisfy 
all requirements, we stop. Otherwise, we lengthen the interval by changing 
B, such that the overlap with the previous row increases. In order to keep the 
total change small, we neglect this approximation and close the row, if there is 
no overlap at all. 

The segmentation resulting from Algorithm |2] satisfies the constraints (i)- 
(iii) and the connectedness, but may contain segments that do not have the 
minimum total field height /. 

4.3.2 Segmentation for ASAS 

The segmentation step for ASAS differs a little bit from the ASS segmentation. 
For ASS, we always find a segment that has its first nonzero row exactly in 
that row where the current matrix B has its first nonzero row. Going through 
the rows, we add further ones to the segment if the current sliding window 
interval overlaps with the previous row. 

Computational tests have shown that for ASAS a different technique makes 
sense because the vertical criterion (iv) plays a role. In detail, whenever we 
decide for Z^+i < li (respectively r^+i > r^), this immediately implies k+k < li 
(respectively r^+fc > r.;) for k = 2, ... ,§2- Additionally, as we need at least gi 
consecutive ones in each row, we also know r^+fc > max(Zi — 1, Z^+i + gi — 1) 
(respectively k+k < min(ri + l,ri+i — gi + 1)). Therefore, we use a matrix 
S = (sij) to store unavoidable ones, i.e. if k+i < U, we put 

Si+fej = 1 for 2 < fc < 32, k+i < j < raa.x{k - 1, k+i + gi - I) (2) 

and if r^+i > r^, we define 

Sj+fcj = 1 for 2 < A; < 32, min(ri + 1, r^+i - gi + I) < j < n+i (3) 

As we require consecutive ones in each row, we also put 

s.,+k,j = 1 for 2 < fc < 32, h+i <j< n+i (4) 

if li+i < U and r,;+i > r,; at the same time. 

We also have to take care that the covered regions have a vertical minimum 
size. If li+i > li (respectively r^+i < ri), we analogously put unavoidable 
zeros into our matrix S using the corresponding rules to ^ and ||3]l. (UJl is not 
necessary here, as zeros do not have to be consecutive in the rows. 

Example 1. Let gi = 32 = 3 and let the (*)-entries of the matrix above 
denote the open bixels for row 1 and 2. Before choosing the leaf positions for 
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row 3, there are some unavoidable ones and zeros that have to be respected. 



/ 



Thus, the choice of the leaf positions in one row produces unavoidable ones 
or zeros in other rows. Our algorithm will choose k and r^ such that the total 
change of this row and the corresponding unavoidable ones is minimal. There- 
fore it might happen, that we do not use the first nonzero row and we also do 
not use the sliding window technique anymore, as the minimum vertical gap 
constraint (iv) prohibits so many leaf positions that we can compare the re- 
maining ones with regard to the resulting total change. Algorithm |4] and |5] in 
the appendix show the corresponding segmentation steps. 

The idea behind this algorithm is, that for each segment and in each row 
we look at all feasible leaf positions. For each pair {I, r), we compute the value 
benchmark, which is the difference between the number of positive entries and 
the number of nonpositive entries in this row as well as in the corresponding 
unavoidable ones. The larger this value is, the better the pair (/, r) suits to the 
segmentation. The unavoidable zeros are not taken into account because it is 
not necessarily bad if an entry hij > is closed, as this entry can be part of 
the following segments. We close a row, if the corresponding optimal value of 
benchmark is zero. 

Finally, our procedure leads to a segmentation satisfying (i)-(iv) and the 
connectedness, but not necessarily satisfies the minimum total field height /. 

Remark 4. One might argue that it is possible for ASS to compute exactly 
the segment, that reduces the total change in this step as much as possible. 
For example, one can consider a layered digraph with m layers of nodes. In 
layer i G [to], we have nodes {i,l,r) representing feasible leaf positions and 
we draw an edge between (i, I, r) and (i + 1, /', r'), if the combination of these 
two leaf positions satisfies all constraints. Furthermore, we draw edges from 
a source to all nodes in the first layer and from all nodes in the last layer to a 
sink. The edge weights are just the total change reductions caused by the leaf 
positions of the end node of the edge and zero for all edges whose end node is 
the sink. For a detailed description of the graph see LSJ. The optimal segment 
can then be found by shortest path computation in the digraph. But indeed, 
such a choice is not a good idea because reducing the total change as much 
as possible leads to badly decomposable residual matrices. For example, for 
gi ^ 2, b — (1, 2, 1) would be reduced by (1, 1, 1) and the residuum (0, 1, 0) is 
badly decomposable. Our used sliding-window-technique is better and leads 
to (1,2,1) = (1,1,0) + (0,1,1). 

For ASAS the constraints are too complex anyway to compute the optimal 
segment in one step, that reduces the total change by a maximum value. 
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Figure 3: The first and the second combination step for ASS with gi 
black rectangle produces an increase of the total change. 



2. The 



4.4 Combination of fields 

Given two segments S = ((^i, ri), . . . , (Z„, r„)) and S" = {{l'i,r[), . . . , {l'^, r'„J), 
let us consider the open regions zi, . . . , 12 and i[,. . . ,i'2of S and 5", precisely 

{ii, . . . , ^2} = {« € H : li<n+ 1} (5) 

{z'i,...,z;,} = {zgH : Z:<n+1} (6) 

If i[ = ^2 + 1 and min(ri2 , r^, ) — max(^i2 ,1'^,) > gi — 1, we merge S and 5" and 
get one new segment S" with 



if Sij = 1 or Sjj — 1, 



4 - '' 
*"' 1 0, otherwise. 

We iterate this step, until no two segments can be merged by this procedure. 
Obviously, this step does not affect the total change. 

Afterwards, we compute a second combination step and merge segments 

if i[ = Z2 + 2 and min(ri2 , r^, ) — max(/j2 ,iy} > gi — I- This means, there is 
only one closed row between the two segments. We compute leaf positions 

/ e [max(Zi2,/^,),min(rj2,r^, ) - gi + 1] and r e [I + gi - l,min(ri2, r^, )] such 

that putting ones to the interval [/, r] in row J2 + 1 produces the smallest total 
change with respect to the current approximation matrix. 

Again, we drop S and 5" out of our segmentation and this time add S" with 




Again, we iterate this procedure, until no more such merges are possible. 
Obviously, this second combination step affects the total change, as we increase 
the approximation matrix by adding ones to the segments. After both of our 
combination steps, all segments still satisfy (i)-(iii), while even more segments 
satisfy (v) now. Another positive consequence is a reduction of the total num- 
ber of (not necessarily pairwise different) segments, called the Delivery Time. 
The combination step is demonstrated in Figure 3. 

For AS AS, we only combine two segments according to one of the two steps 
described above, if the criterion with the minimum vertical zeros (iv) is not 
violated after the combination. 

4.5 Make-two-of-two 

We now define a substitution step S + S' = S + S'. For AS AS, the same step is 
computed if S and S' still satisfy (iv). 
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Figure 4: The Make-two-of-two-step for ASS with 51 = 2 and f = 4:. 

Let us again consider a segment S* = ((^1, ?'i), ■ • ■ , ('m, J'm)) and its open 
region ii, . . . , ^2 defined by (|5). If 5 violates (v), we call S a critical segment. 
Now we check whether we find a segment S" = {{l[, r[), . . . , (/^, r^)) with 
open region i[, . . . ^i^ defined by (|6) such that 



«2 



i'l > / - 1 and i'a - n > / - 1. 



Note that in this case, we have i'l < ii and 12 > 12 due to 12 — ii < / 
the set of open rows of 5 is a subset of the set of open rows of S". If 



miii(rij,r^^_i) - niax(/j^ 



we substitute S and S' by segments 5 = ((^i, fi), . 
S' ^ (([;, fi), . . . , (4, f;„)) defined as follows 



/■,-i)>5i-l, 

5 V^rn? ^?7i 



1. Thus, 
(7) 



))and 



iii,ft) 




if i < ii, 
if i > ii, 



(^^,r: 






(n + 1,0) 



if i < ii, 
Hi > ii. 



The result is that we add the upper part of segment S' to segment S in order 
to enlarge S, while 5' remains sufficiently large. If is not satisfied, we can 
have a second try and check whether 



niin(ri. 



_i) -max(/,, -C+i) >.gi -1 



(8) 



If this condition is true, we can analogously add the lower part of S' to S and 
close all rows > 12 + 1 of 5'. The Make-two-of- two-procedure is illustrated in 
Figure 4. 

The make-two-of-two-step is computed for all critical segments. If we find 
a partner, we compute the substitution immediately. If no partner is found 
throughout the segmentation, the segment S is dropped and stored in a new 
list of critical segments violating (v). Whereas the substitution with a partner 
does not affect the total change, the elimination of segments from the segmen- 
tation that find no partner leads to an increase of the total change. 



4.6 Handle Critical Segments 

Let S = ((/i, ri), . . . , {Im, rra)) be a critical segment stored in the make-two-of- 
two-step and let ii , . . . , 12 be its open region defined by (|5). We try to combine S 
with a partner S' = {{l'i,r[), . . . , (Z^, r'„J) from the segmentation with open re- 
gion i[, . . . ,{'2 defined by (|6j if this combination step causes lower total change 
then simply omit S. 

More precisely, if ij > 12 and i[ < ii, we define a new segment S" = 
(a'/, <),..., (CO) with 



/ill II \ 
ih > r, ) 



(min(?i,/9,max(rj,r^)), iiii <i <i2, 
(/^,r9, otherwise. 
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Figure 5: The Handle-critical-segments-step for ASS with gi — 2 and / = 4. 
The black rectangle shows the combining bixels increasing the total change. 

The resulting segment S" emanates from S' by attaching certain (or all) 
parts of S and possibly also adding connecting ones. We now drop 5" from 
our segmentation and S from the list of critical segments and add S" to the 
segmentation, if the total change that is caused by this decision is smaller than 
simply omitting S. If no partner for 5* is found, we simply delete S from the 
list of critical segments and accept the corresponding total change. 

For AS AS, we execute the step if (iv) is not violated. 

4.7 Total Change Improvement 

As the previous segmentation steps are indeed not optimal, it may happen, that 
some leaf positions can be increased (respectively decreased) without violating 
the constraints and the total change becomes smaller. 

Let B be the approximation matrix that corresponds to the current segmen- 
tation. The improvement step is computed using Algorithm |6] which is de- 
scribed in the appendix and adjusts the segmentation such that B gets closer 
to A with regard to the total change. 

5 Clinical case 

A clinical case was set-up to examine the efficiency of our proposed segmen- 
tation algorithms. For a patient with cancer of the right breast, electron irra- 
diation plans using various segmentation settings and different optimization 
constraints were created with a self-designed IMRT optimization programme 
based on our previous studies |9l- The planning target volume was the right 
breast, which should receive a total dose of 50.4 Gy (1.8 Gy per fraction). In 
addition, the target volume should be covered by the 95% isodose line (95% of 
the prescribed dose). The ipsilateral lung was considered to be organ at risk. 

The optimization programme provides simultaneous optimization of beam 
orientation, energy and intensity for dose delivery with an add-on MLC for 
electrons (Euromechanics, Schwarzenbruck, Germany) presented in Figure [6] 
and IIT2IIT3I . Electron dose calculation was performed by Monte Carlo simula- 
tions with the treatment planning system Pinnacle from Philips (Version 8.1s). 
Final dose calculation of the treatment plans was conducted using a dose grid 
size of 3 mm and a dose calculation uncertainty of 2%. The following optimiza- 
tion steps are necessary to generate an electron IMRT plan: 

1 . Simultaneous optimization of beam orientation, energy and intensity: A 
set of radiation incidence angles (typically 10-15) is determined given by 
table and gantry angles |9, 24|. For each configuration, the algorithm cal- 
culates the optimal fluence distribution, given by a nonnegative integer 
matrix. 
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2. The intensity matrices are approximately decomposed into a superposi- 
tion of allowed segments such that the deviation between desired and 
actual fluence is minimal. The result is a set of segments, each of them is 
given by the corresponding MLC leaf positions and its dose weight. 

3. The segments from step 2 are treated as candidates for the treatment plan. 
In a third optimization step, the dose of the candidates are calculated for 
all beam energies and then optimized for a given weight proportion be- 
tween best target coverage and minimum dose to critical organs in order 
to find the final set of segments with optimal beam energies and their 
corresponding monitor units. 

In this paper, we have focused on step 2 and introduced optimization algo- 
rithms for an approximate decomposition of intensity matrices. 

Until now, the segmentation step consisted in exactly decomposing the in- 
tensity matrices using all deliverable segments. In our approach, we admit 
a decrease in the decomposition accuracy in order to obtain segments which 
satisfy the dosimetric and technical constraints. Step 3 justifies the approxima- 
tion approach in Step 2, as a larger approximation error does not necessarily 
result in a suboptimal treatment plan. Indeed, larger segments produce homo- 
geneous dose distributions and thus, the same final fluence can be generated 
using fewer larger segments. The acceptability of a treatment plan is decided 
after step 3 by means of dose volume histograms (see Section [SJl and a plan is 
only presumed if the required dose constraints are not exceeded. Therefore, 
the danger of cumulative deviation in the approximation step does not really 
exist, as the computed segments are just candidates for the treatment plan that 
pass through a further optimization step. 

6 Results 

At first, we compare electron IMRT plans created with different segmentation 
settings for the clinical case prescribed in Section|5] The comparison was con- 
ducted using two different optimization settings: one setting to achieve with 
a better dose coverage of the breast and the other one to reach a better spar- 
ing of the lung. Note, that the segmentation settings refer to parameters in the 
segmentation step, which is discussed here, whereas the optimization settings 
play a role in the final optimization step that is not part of this paper. Finally, 
we give a detailed evaluation for the results of the decomposition step. 

A treatment plan with a segmentation setting xyz uses the decomposition 
algorithm with a minimum total field height oi f — x, a minimum separation 
constraint and row overlap oi gi = y and a minimum vertical gap of 92 = z. 
The decomposed matrices vary in their vertical size m and their horizontal 
size n, as they describe only parts of the beam head where the target volume 
is located. Thus, in practice, the overtravel parameters bi and br will depend 
on the positioning of the matrix and are put individually for each matrix. Our 
electron MLC is capable of shifting the leaves edges to 3/4 of the radiation field. 

The plan quality was evaluated by means of dose volume histograms that 
indicate the amount of dose delivered to a certain volume of the patient (here: 
the right breast and the right lung). Thus, dose homogeneity in the target vol- 
ume and dose exposure to the organs at risk can be examined. In Figure [T) the 
dose volume histograms for both optimization settings demonstrate that al- 
most identical dose distributions can be achieved using smaller or larger mini- 
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Figure 6: Left: Add-on MLC for electrons mounted on a conventional Siemens 
accelerator. Right: Dose distribution of an electron IMRT plan consisting of 
26 MLC fields achieved through segmentation setting 442. The corresponding 
dose volume histogram is shown in Figure [7] (left). The setting 442 is given by 
a minimum total field height of 4, a minimum separation constraint and row 
overlap of 4 and a minimum vertical gap of 2. 



mum MLC openings (cf. setting 111 and 441). In fact, the treatment plan could 
be slightly improved by use of a minimum vertical gap parameter of 2 which 
avoids single leaf openings and closings. 

Table[T]and|2]illustrate the main benefit of our approach, as identical results 
can be achieved with approximately two thirds fewer segments and a signifi- 
cantly smaller number of monitor units by use of dosimetric constraints. As a 
result, the dose delivery is more efficient and less time consuming. The min- 
imum number of segments is reached for setting 442 and computational tests 
have shown, that this configuration produces the optimal results. As the leaf 
width is 0.7 cm, fields with a horizontal and vertical height of 4 bixels have a 
size of approximately 3 cm x 3 cm and this confirms our dosimetric constraint 
of 3 cm X 3 cm minimum segment size (cf. Figure [l). It can be also demon- 
strated that minimum segment sizes greater than setting 442 do not necessarily 
result in fewer segments (cf. Table [1), although the number of segments in Ta- 
ble|2]is slightly lower for setting 552. For both optimization settings, the dose 
volume histograms were considerably better when using minimum segment 
sizes smaller than setting 552. 

It is important to underline that the number of segments and the number 
of monitor units in Table[T]and|2]belong to the final IMRT plan and result from 
the third optimization step and not from the decomposition step of our algo- 
rithms. In fact, the monitor units have another scale here and are not directly 
comparable with the delivery time from the segmentation. In contrast, the to- 
tal change information stems from the decomposition step. Note, that the total 
change of the segmentation itself is not a significant quantity, because if the 
matrix entries are large, a larger total change is acceptable. Therefore, we com- 
pute the total sum of entries for each intensity matrix and then calculate the 
relative total change which is the ratio between total change and total sum of 
entries. The smaller the relative total change, the better is the decomposition. 

For the detailed evaluation of our algorithms, we use a set of 264 clinical 
intensity matrices that originate from electron treatment plans for different pa- 
tients and beam angles. The matrices are produced during the optimization 
step 1 of the treatment planning that was introduced in Section |5] and uses the 
algorithm from |9|. Exemplarily, we compute segmentations for the settings 
/ = 3, gi = 3 and 92 = 1 as well as / = 4 and gi = 92 = 2. The values of 
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Figure 7: Dose volume histograms for settings 111, 441 and 442 to demonstrate 
that identical results could be achieved when using greater minimum segment 
sizes (cf. Ill and 441) and segment shapes without single vertical leaf openings 
and leaf closings (cf . 441 and 442). The setting xyz is given by a minimum total 
field height of x, a minimum separation constraint and row overlap of y and 
a minimum vertical gap of z. Left: Dose volume histogram for optimization 
setting to reach a better dose coverage of the breast volume. The resulting 
dose distribution for setting 442 is presented in Figure [6l Right: Dose volume 
histogram for optimization setting to provide a better sparing of the lung. 



the overtravel parameters are also produced in the pre-segmentation step. The 
results are shown in Table |3] and |4] and demonstrate how much total change 
is caused respectively avoided by the steps of the algorithms. The overtravel- 
approximation and MSC-approximation lead to a certain total change of the 
matrix that is put into the segmentation step. As an exact decomposition in the 
segmentation step is impossible, the total change increases here again. Both 
the combination step and the make-two-of-two step try to eliminate segments 
not satisfying the parameter / and again cause some total change. Finally, the 
last two steps of the algorithm improve the performance and reduce the ap- 
proximation error as much as possible. One can see that the combination step 
and the make-two-of-two-step are performed more often for ASS, as for ASAS 
the vertical gap ensures that the fields already have a reasonable size after the 
segmentation step. Of course, the larger the parameters and thus the minimum 
field size, the larger becomes the total change. 

The first column in Table |3] and |4] gives the average results over the 264 ma- 
trices, while the second (respectively third) column represent the single results 
for the matrix with the smallest (respectively largest) relative total change. Ho- 
mogeneous matrices with large nonzero areas can be decomposed quite well, 
while matrices with only few nonzero entries that do not span connected areas 
lead to unacceptable results. As a treatment plan is a superposition of several 
intensity profiles from different beam angles, the approximation errors balance 
each other and lead to applicable treatment plans as described above. Fur- 
thermore, a certain part of the total change is unavoidable if one requires the 
constraints (i)-(v), e.g. the total change after MSC-approximation is a good 
lower bound for the achievable total change. All in all, taking the vertical gap 
52 into account increases the total change while reducing the number of used 
segments and the monitor units. 
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0.04 


221 


0.16 


222 


0.34 


331 


0.22 


332 


0.36 


333 


0.40 


441 


0.30 


442 


0.41 


443 


0.45 


444 


0.49 


551 


0.35 


552 


0.46 



Table 1: Segmentation results of IMRT plans using different decomposition 
settings for optimization setting to achieve the best target coverage. Setting 
xyz means f = x (minimum total field height), gi = y (minimum separation 
constraint and row overlap) and 52 = ^ (minimum vertical gap). 

Setting Mean Relative Total Change Number of Segments Number of Monitor Units 

90 85586 

79 65677 

54 63461 
49 38598 
32 21265 
40 36789 

55 34262 

26 12860 
28 13119 

27 16899 
40 19162 
30 10337 



Table 2: Segmentation results of IMRT plans using different decomposition set- 
tings for optimization setting to reach less dose to the lung. Setting xyz means 
f ~ X (minimum total field height), gi = y (minimum separation constraint 
and row overlap) and 92 = z (minimum vertical gap). 

Setting Mean Relative Total Change Number of Segments Number of Monitor Units 

78 59399 

71 46627 

51 42411 

52 25896 
29 16700 

36 24957 
45 23509 
28 10200 
28 9995 
27 11318 

37 13787 
25 8832 



7 Conclusion 

In the present study, dosimetric and technical constraints have been taken into 
consideration in intensity-modulated radiation therapy (IMRT). A set of 5 pa- 
rameters has been introduced, two of them for modelling the leaf overtravel 
constraint, the other three to ensure a minimum field size and to avoid thin 
field shapes. We proposed algorithms for approximate segmentation of inten- 
sity matrices using segments that satisfy the constraints. We basically distin- 
guish between two approximation problems depending on whether the ver- 
tical gap parameter is considered or not. The objective function of the opti- 
mization is the deviation between the desired and the approximated intensity 
profile that has to be minimized. The segmentation step is part of an IMRT op- 
timization process which was examined by comparisons of dose volume his- 
tograms of treatment plans with small and large segments as well as with and 
without thin segment shapes. The histograms show that the use of larger seg- 
ments results in equal IMRT plans with fewer segments and monitor units re- 
spectively. Although the approximation error of the segmentations rises with 
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0.16 


222 
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331 


0.22 


332 


0.36 


333 


0.40 


441 


0.30 


442 


0.41 


443 


0.45 


444 


0.49 


551 


0.35 


552 


0.46 



Table 3: Numerical results for ASS with parameters (?i = 3 and / = 3. 





Average 


Rel. TC min. 


Rel. TC max. 


m 


17.55 


20 


11 


n 


21.44 


19 


13 


Total sum of entries 


886.4 


1329 


152 


Total change 


82.24 


12 


64 


Delivery time 


15.39 


13 


3 


Number of segments 


15 


13 


3 


TC after Overtravel-Approximation 


26.34 


6 


75 


TC after MSC- Approximation 


57.09 


11 


88 


TC after Segmentation 


74.33 


13 


91 


TC change after combination 


80.96 


28 


91 


TC after make-two-of-two 


93.77 


28 


101 


TC after handle critical segments 


88.72 


28 


94 


TC after improvement 


82.24 


12 


94 


Combinations 


9.19 


5 





Successful make-two-of-two 


2.78 








Relative total change 


0.16 


0.009 


0.62 



Table 4: Numerical results for ASAS with parameters gi = g2 — 2 and f = 4:. 





Average 


Rel. TC min. 


Rel. TC max. 


m 


17.55 


23 


22 


n 


21.44 


26 


28 


Total sum of entries 


886.4 


997 


384 


Total change 


157.7 


46 


269 


Delivery time 


7.69 


4 


4 


Number of segments 


6.82 


3 


4 


TC after Overtravel-Approximation 


26.34 


7 


89 


TC after MSC- Approximation 


44.39 


7 


134 


TC after Segmentation 


144.8 


46 


267 


TC change after combination 


144.8 


46 


267 


TC after make-two-of-two 


159.2 


46 


270 


TC after handle critical segments 


158.8 


46 


270 


TC after improvement 


157.7 


46 


269 


Combinations 


0.24 








Successful make-two-of-two 


0.57 








Relative total change 


0.27 


0.05 


0.70 



increasing minimum field size, equivalent or even better dose distributions 
could be achieved. Concluding, this first approach to approximated segmen- 
tation in IMRT planning shows the potential of these ideas and there is a need 
for further research in related approximation problems. 
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8 Appendix 



Algorithm 1 LOC-left 



iput: vector t; 




for i — min to max do 




tcii = \vi - i\ 




end for 




tcij — oo for all i G [min, 


max] , j 


for J = 2 to fc do 




for i = min to max do 




for i' = min to i do 




if icij_i + i^j — i' 


1 < icj/j 


tCi'j — tCij^i + 


h-^'\ 


prcifj = I 




end if 




end for 




end for 




end for 





]>l 



then 



opt Ulll^ie [min, maa;] ''Ci,fc 

Let Zopi be one of the indices with tci^^^,k ~ opt. 
for j ~ k down to 1 do 

If still j > 1, then iopt = prci^ptj 
end for 
Output: vector w 



Algorithm 2 Segmentation ASS 



Input: matrix B 
while B 7^ do 

while hi = do 

li = n + 1; Vi = 0; i = i + 1; 
end while 

Choose k and r^ for row i such that Vi — li> gi — 1 (see Remark 2). 
while Ti — li > gi — 1 and bi-|_i 7^ do 

Find interval ASS: Choose Z^+i and r^+i for row « + 1 such that the 

overlap with the previous row is > gi or close row i + 1. 

i = i + 1; 
end while 

Close all remaining rows with k = n + 1 and r^ = 0. 
for I = 1 to TO do 

for j — li to Ti do 

end for 
end for 

Store ((^1, ri), . . . , (/,„, Tm)) in the segmentation. 
end while 
Output: Segmentation 
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Algorithm 3 Find interval ASS 



Input: i,li,ri 
finish = 0; 
while finish = do 

start = min{j | &i+i j- > 0}; 

end = minlj | bi+ij > 6i+ij+i}; 

I — max(siart, li); 

r — min(en(i, r^); 

ii r — I > gi — 1 then 

/i+i — start; r^+i — end; finish = 1; 
else if r Z < then 

/i+i ~ n + 1; r^+i = 0; finish — 1; 
else 

if start > li then 

&i+l,start-l = bi+i^start-1 + 1; 

else 

bi+l^end+l = bi^i^end+1 + 1; 

end if 
end if 
end while 
Output: h+i, ri+i 



Remark 2. The interval [I, r] for the first open row is computed analogously 
to the interval for the other rows, only ignoring the overlap constraint and 
instead requiring r — I > gi — 1. 
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Algorithm 4 Segmentation ASAS 



Input: matrix B 
while B y^Odo 

i = 1; Sij = for all i and j; 
while 6j = do 

li ^ n + 1; Vi = 0; i = i + 1; 
end while 
start = i; closed = 0; 
if start > n — 92 + 1 then 

closed = 1; 
else 

Choose li and r^ for row i such that ri — k > gi — 1 or close row start 

(see Remark 3). 

if row start is closed then 

closed = 1; 
else 
while r.i — It > gi — I and 5i_|_i 7^ do 

Find interval ASAS: Choose k+i and r^+i for row i + 1 such that the 
overlap with the previous row is larger than gi or close row i + 1. 
i = i + 1; 
end while 

for all remaining rows i do 
if Sj ^ then 

li =index of the first one in Sj; r^ =index of the last one in Sj; 
else 

li = n + 1, ri = 0; 
end if 
end for 
end if 
end if 
if closed — 1 then 

bstart = 0; 

else 

for i = 1 to m do 
for j = li to ri do 
if bij > then 

end if 
end for 
end for 

Store ((/i, ri), . . . , {Im, rm)) in the segmentation. 
end if 
end while 
Output: Segmentation 
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Algorithm 5 Find interval ASAS 



Input: i + l,li,ri, current matrix S 
if Si = theri 

vi = n + 1; V2 = 0; 
else 

fi = index of the first one in si; V2 = index of the last one in Sj; 
end if 

ii = index of the first positive entry in bj; 
^2 = index of the last positive entry in 6^; 
opt == 0; k+i = n + 1; r^+i = 0; 
for / — min(wi,ii) to m.m(bi,vi) do 

for r ~ niax(Z + gi — 1,V2, br) to niax(u2, ^2) do 

if niin(r, r^) — niax(Z, k) > gi — 1 and not (l < li and « + 1 > m — 52 + 1) 
and not (r > r^ and i + 1 > to — 92 + !) and no unavoidable zero is 
between I and r then 
benchmark = 0; 

for all entries {k,j) of unavoidable ones with « + 2 < fc < i + 92 do 
if bkj > then 

benchmark — benchmark — 1; 
else 

benchmark = benchmark + 1; 
end if 
end for 

for j — I tor do 
if bi+ij > then 

benchmark = benchmark — 1; 
else 

benchmark ~ benchmark + 1; 
end if 
end for 
if benchmark < opt then 

opt ~ benchmark; li^i = /; r^+i = r; 
end if 
end if 
end for 
end for 

Put tinavoidable ones or zeros corresponding to the decision for Z^+i and r^+i 
into S. 
Output: k+i, r^+i, current matrix S 



Remark 3. The choice of {I, r) in row start is again computed analogously, 
only ignoring the overlap constraint with the previous row. 
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Algorithm 6 Total Change Improvement 



Input: Segmentation 

for all segments do 

for i = 1 to TO do 

i{li> 1 and bi^i.^i < ai^i^^i 

(for ASAS: and decreasing k will not violate (iv)) then 

H H J-; 

end if 

if ri < 1 and bi^n+i < (^i,ri+i 

(for ASAS: and increasing r^ will not violate (iv)) then 
ri=ri + 1; 
end if 
if bi^i- > Qi^i^ and increasing li will not violate (iii) 

(for ASAS: and increasing k will not violate (iv)) then 
H — H ' r, 
end if 

if bi,ri > ai,ri arid decreasing r^ will not violate (iii) 
(for ASAS: and decreasing r^ will not violate (iv)) then 
ri = ri- 1; 
end if 
end for 

Update approximation matrix 
end for 
Output: Segmentation 
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